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Abstract 

, We use the recently developed Center for Radiative Shock Hydrodynamics (CRASH) code to numerically simulate 

' laser-driven radiative shock experiments. These shocks are launched by an ablated beryllium disk and are driven down 
xenon-filled plastic tubes. The simulations are initialized by the two-dimensional version of the Lagrangian Hyades 
D ' code which is used to evaluate the laser energy deposition during the first 1 . 1 ns. The later times are calculated with the 
, CRASH code. This code solves for the multi-material hydrodynamics with separate electron and ion temperatures on 
QQ ' an Eulerian block-adaptive-mesh and includes a multi-group flux-limited radiation diffusion and electron thermal heat 
conduction. The goal of the present paper is to demonstrate the capability to simulate radiative shocks of essentially 
three-dimensional experimental configurations, such as circular and elliptical nozzles. We show that the compound 
shock structure of the primary and wall shock is captured and verify that the shock properties are consistent with 
order-of-magnitude estimates. The produced synthetic radiographs can be used for comparison with future nozzle 
experiments at high-energy-density laser facilities. 
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In the experiments of the CRASH project, high-energy-density plasma flow is driven by a beryllium disk irradiated 
^ ' by the Omega laser beams. The ablation of part of the beryllium results in a rocket launch of the remaining beryllium 
, driving a strong shock-wave through a xenon-filled plastic tube. The shock heats the xenon gas to sufficiently high 
temperature so that it will ionize and create free electrons. Directly behind the shock the electrons and ions will 
equilibrate due to the Coulomb collisions. This will heat the electrons. The shock is fast enough so that the electrons 
have to emit radiation behind the shock in order to satisfy the energy balance equation. This results in a radiative 
, cooling layer 111]. The transport of the radiation ahead of the shock will heat and ionize the xenon in the radiative 
precursor. Some fraction of the radiation in the upstream region will also transport sideways and strike the plastic 
tube. These photons can heat the plastic ahead of the primary shock. The ablated plastic that moves then inward 
compresses the xenon resulting in a so-called wall shock QUI]. It is the interplay between the primary and wall shock 
that is of interest for the CRASH project. 
■ To model the laser-plasma physics and laser energy deposition in the radiative shock tube experiments, we use 

^ , H2D, the 2D version of the Hyades code iQ], which contains a built-in laser package. H2D is a Lagrangian radiation- 
hydrodynamics code that utilizes the axisymmetry. Hyades is capable of tracing rays in 3D; the runs shown below 
used 2D ray tracing for the laser energy deposition. We use the H2D code to simulate our experiments for the first 
1.1 ns (sometimes up to 1.3 ns), the time of the full width half maximum (FWHM) 1 ns laser pulse including ramp- 
up and ramp-down time. The spatial profile of the beam is determined using a model for the irradiance pattern of 
the laser beams in typical experiments. At 1.1 ns the experiment is in a regime that is well described by radiation- 
hydrodynamics so that the simulation can be continued with the CRASH code jstl instead. This code solves for the 
multi-material hydrodynamic equations with a multi-group flux-limited diffusion model for the radiation and uses the 
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recently developed Block Adapative Tree Library (BATL) We are currently constructing a laser package in our 
simulation code as an alternative to H2D and the progress in that work will be reported elsewhere. 

The main goal of the CRASH project is to assess and to improve the predictive capability of a simulation code, 
based on the combination of experiments, simulation studies and statistical analysis. The baseline experiment is a 
laser-driven shock launched through a xenon-filled straight tube. We plan to perform an experiment in which the tube 
geometry is changed to contain a circular nozzle. These experiments and the comparison with the simulations will 
demonstrate how well the numerical code can predict the shock properties of the following experiment in which a 
wide circular tube tapers into an elliptical tube. In this project we plan to analyze the radiative shock structures in 
nozzle shock tubes in the context of predictive capability and uncertainty quantification of the simulations [1. The 
demonstration that such numerical simulations are feasible is the topic of the present paper 

This paper is organized as follows. Section |2] describes the radiative shocks produced in nozzles with circular 
cross-section. We contrast the results with order-of-magnitude estimates based on physical arguments. Section |3] 
demonstrates the radiative shocks in nozzles with an elliptical neck. This simulation is fully three-dimensional. The 
paper concludes with Section]?] 



2. Circular nozzle 

A laser pulse irradiates a 20 /im thick beryllium disk with 0.35 fim wavelength light for the FWHM duration of 1 ns 
and with a laser energy deposition of 4kJ. For the laser spot size we use a FWHM 800 jum diameter which is smaller 
than the 1200 ;um diameter of the tube. We are currently constructing a laser package in the CRASH code, but the runs 
reported here use the Lagrangian radiation-hydrodynamics code Hyades 2D (H2D) |4] to evaluate the laser energy 
deposition during the first 1.1 ns which includes the laser ramp-up and ramp-down time. While our aim is to simulate 
radiative shocks in nozzles, the first 1.1 ns is however simulated in a straight tube for convenience. Calculating the 
first 1.1 ns with a straight tube and then transforming this tube into a nozzle is physically justified as long as the taper 
and shaft of the nozzle does only alter the radiative precursor The justification originates from the observation that the 
radiative transport from the precursor back through the shock front is negligible [8]. The considered straight tube does 
have a cylindrical polyimide wall of 100 fim thickness and inner radius of 600 yum filled with xenon with mass density 
p = 0.0065 g/crn^, while there is vacuum outside. The beryllium disk is immediately to the left of jc = 0, where x is 
the coordinate along the tube. The laser light will come in from the negative x direction. We place a gold washer next 
to the beryllium disk to protect the outside of the plastic tube from the laser light and use acrylic in between the gold 
and polyimide tube. 

After 1.1 ns of simulation time the output of H2D is used to initialize the Eulerian CRASH code as described in 
Appendix A The straight tube of H2D is transformed into a nozzle, see the left panel of Fig. lT]for one quarter of 



the nozzle domain. This nozzle changes cross-section in the following way. For x < 500 i^m we do not modify the 
tube as defined in H2D. For x > 750//m we shrink the tube diameter from 1200 //m to 600 ;um and the polyimide 
wall thickness corresponding reduces to 50 yum. Between 500 jim and 750 jim the tube diameter and wall thickness are 



linearly shrinking (In the notation of Appendix A we use xq = 500 //m, xi - 750;um, Sy -1/2 and - 1/2). The 
vacuum outside the nozzle is replaced with low density polyimide to avoid the otherwise zero mass density. The left 
panel of Fig. JT] shows the materials that are present in the simulation in color: beryllium (blue), polyimide (green), 
acrylic (red), and gold (yellow). The xenon inside the nozzle is for convenience not shown in this figure. 

The computational domain size is -150 < x < 3900, < y < 900 and < z < 900 in microns, where y and z are 
the two directions transverse to the nozzle. The simulation is performed with an effective resolution of 2560x5 12x5 12 
grid cells using two levels of refinement. The effective cell sizes are therefore approximately 1 .6 fim along the tube 
and 1 .8 ;um in the two transverse directions. The domain is decomposed in 4 x 4 x 4 grid blocks. The mesh is refined 
at all interfaces that involve xenon or gold. To capture the shock front and the cooling layer, the mesh is also refined 
where the xenon density exceeds 0.02 g/cm^. Grid blocks are also refined if these criteria are satisfied in the ghost 
cells. 

We perform the simulations with the CRASH radiation hydrodynamics code [5]. This code solves for the radiation 
hydrodynamic equations in three operator splitting steps: (1) an explicit time step of the hydrodynamic equations 
using a shock-capturing solver, (2) a linear advection of the radiation bins in frequency-logarithm space, and (3) an 
implicit solve of the radiation diffusion, heat conduction, and energy exchanges. This code is continuously undergoing 



2 



Figure 1: The materials in the 3D circular nozzle experiment at 1.1 ns (left panel) and 13ns (right panels). Only one quarter of the shock tube is 
shown. The materials are beryllium (blue), polyimide (green), acrylic (red) and gold (yellow). The xenon inside the tube is not colored, but the 
border of the xenon volume is indicated by the white surface. The primary shock at 13 ns is depicted by a black iso-surf'ace. 



improvements. One such improvement is the resolution change treatment for the radiation diffusion and electron heat 



the original scheme in [SJ was used instead. We use for the hydrodynamic part of the equations the HLLE scheme 
with a Courant-Friedrichs-Lewy number of 0.8 and the generalized Koren limiter with y6 = 3/2. For the radiation, we 
use the multi-group flux-limited diffusion model with 30 groups. The photon energy range is 0.1 eV to 20keV that is 
logarithmically distributed over the groups. The radiation diffusion, heat conduction and energy exchanges are solved 
with the split (decoupled) implicit solver of [5] using the conjugate gradient method with a Schwarz-type Incomplete 
Upper-Lower (ILU) preconditioned 

Due to the symmetry in the problem we only need to simulate one quarter of the nozzle domain (y > and z > 0) 
and use reflective boundary conditions at y = and z = 0. For all other boundaries of the domain we use extrapolation 
with zero gradient. For the radiation we use zero albedo boundary conditions. 

We simulated the shock evolution from 1.1 ns to 13 ns physical time. It took a little over three days to compute on 
1000 cores of the HERA supercomputer at the Lawrence Livermore National Laboratory. The number of cells in the 
computational domain increased from 26.5 million at the beginning to about 38 million near the end. The 3D material 
identity at 13 ns is shown in the right panel of Fig. [1] The beryllium has moved into the nozzle and is like a piston 
driving a shock in the xenon. This shock is indicated by a black iso-surface (of high ion temperature values). The 
xenon itself is not shown but the edge of the xenon is emphasized by a white color to make the xenon entrainment 
between the beryllium and polyimide more clear. We can also see the inward moving polyimide which will lead 
to a wall shock. In the following we will analyze the shock structure in more detail and check if the results are in 
agreement with back-of-the-envelope estimates presented in 

The shock structure in the xy-plane at time 13 ns is shown in Fig. |2] The top left panel is for the material 
identification of beryllium (blue), xenon (black), polyimide (green), gold (yellow) and acrylic (red). The nozzle with 
inner radius of 600 fj.m, taper, and shaft with inner radius of 300 //m are visible. The black lines indicate the resolution 
changes between the grid refinement levels. The top right panel shows the mass density. Part of the polyimide is of 
very low density and represents vacuum. The dense polyimide tube thickness ranges from 50 to 100 fim. The xenon is 
compressed by the beryllium piston flow resulting in a primary shock that is located at x x 1700 fim. For convenience 
we have indicated with black lines where the material interfaces are. 

To check the obtained properties of the primary shock we first determine the beryllium piston velocity. For a 1 ns 
laser pulse of 4 kJ energy, the irradiance on a 800 //m diameter spot size is 8 x lO''* W/cm^. Most of this light will be 




During the simulations of the nozzles these improvements were not yet present and 
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Figure 2: The radiative shock structure in the 3D circulai' nozzle simulation at 13 ns. The plots show in color contour the plasma and radiation state 
indicated in the plot titles as a function of the x and y positions in microns. The colors in the top left panel indicate beryllium (blue), xenon (black), 
polyimide (green), gold (yellow) and acryllic (red). The black lines in this panel show resolution changes, while in the top right panel the lines 
indicate the material interfaces. 
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Figure 3: Zoom-in of the ion temperature near the shock with the .v and y coordinates in microns. 



absorbed in the beryllium, so that the absorbed energy per unit of area during this 1 ns is 8 x 10^ J/cm^. About 20% of 
the beryllium mass, corresponding to 4 yum of the 20 fim, will be ablated by the laser [8]. The ablation efficiency, which 
is the ratio of the kinetic energy of the remaining 16 fim beryllium to the total kinetic energy and exhaust, is to lowest 
order in the ablation percentage equal to the fraction of the mass that initially has been ablated, i.e. the efficiency is 
20%. In reality, however, about half of the incident laser light reaches the beryllium above the absorption region where 
ablation occurs, while the other half is absorbed below the critical density. The actual efficiency from converting laser 
energy to kinetic energy of the remaining 16 //m of beryllium is therefore roughly 10%. The areal mass density of 
16^(m beryllium at 1.8 g/cm^ is approximately 3 x lO^-'g/cm^. The areal kinetic energy density 0.1 x 8 x lO^J/cm^ 
corresponds therefore to an initial velocity of the beryllium of a little more than v = 200km/s. This beryllium will 
launch a shock through the xenon-filled tube. This high velocity is only achieved at early times. The numerically 
obtained shock velocity gradually reduces to a value between 1 lOkm/s and 120km/s at time 13 ns as shown in the X 
velocity plot of Fig. |2] 

The unshocked xenon gas pressure is initially about 1.1 atm while the density is p = 6.5 x lO^^* g/cnP. With 
the above mentioned shock velocity of v a; 1 10 km/s we obtain a xenon post-shock pressure at 13 ns of the order of 
pv^ X 80 GPa. The post-shock pressure at x x 1700 obtained in the bottom right panel of Fig. |2]is in agreement with 
this estimate. The shock wave heats the ions. The ion temperature in the postshock region of a strong shock wave 
with compression ratio equal to k - (y + l)/(y - 1) is for our application approximately H: 



kBTj - AnipV 



jl-l/^ 



(1) 



where kg, nip and A are the Boltzmann constant, the proton mass and atomic mass, respectively. For xenon with atomic 
mass A - 131, a shock velocity of 1 10 km/s, and an adiabatic index of y = 5/3 for ions gives a postshock temperature 
of r, a; 3keV. We find with the numerical simulation an ion temperature of 1.2 keV in Fig. [3] that shows the region 
around the shock. With an effective resolution of 1 .6 jum we do not fully resolve this narrow ion temperature peak, 
resulting in a lower than expected maximum temperature. We find two other temperature peaks as well. One peak is 
near x - 1600 //m behind the shock and another is near y - 250 jim. They overlap with the beryllium-xenon and the 
polyimide-xenon material interfaces, respectively. These spikes are probably the result of a low-order convergence 
rate in the interface treatment. 

For a high density plasma the collision frequency between electrons and ions is large, so that their temperature will 
equilibrate. At the shock, however, the ion temperature jumps so that the ions and electrons are out of equilibrium. 
The ions will heat the electrons via Coulomb collisions and form an equilibration zone directly behind the shock. The 
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electron heating will also increase the ionization in the equilibration zone. Fig. |4]shows the average ionization in the 
xy-plane. We find from the code that the ionization is elevated to about (Z) - 17 at the electron temperature peak in 
the equilibration zone. At this peak the energy of each ion is shared with (Z) electrons. The equilibration temperature 
Teq can be approximated as the postshock temperature for a strong shock under the assumption that the electrons and 
ions are in temperature equilibrium: 

(Z) +1 K 

For representative values of the polytropic index between y - 1.2 and 1 .3 for single-temperature xenon, the estimated 
equilibration temperature is between Teq - 76 eV and 104eV. In reality this temperature will be somewhat lower 
due to radiative cooling. The numerically obtained value is ~ 73 eV in Fig. |2] This electron temperature is not 
the final state in the postshock region since crT^ « 2.8 x lO'^ W/cm^, where cr is the Stefan-Boltzmann constant. 
The incoming kinetic energy flux is however ^ 4.3 x lO" W/cm- based on the velocity at time t = 13 ns 

of V = llOkm/s and xenon density p - 0.0065 g/cm^. There is therefore not enough incoming energy. There 
must be a cooling layer ^ through which the electron temperature falls to the final temperature Tf estimated by 
2o-Tj = * 4.3 X 10" W/cm2. The factor 2 is bee ause the radiative cooling layer emits in both directions 

equally. The final temperature is thus Tf ^ 38 eV. 

The heated electrons are the main energy source for radiation. In Fig. |2]we show the radiation temperature as a 
measure for the total radiation energy density. The photons travel upstream of the shock where they preheat and ionize 
the unshocked xenon in the precursor as depicted in the electron temperature panel in Fig. |2]and Fig. |4] The radiation 
transport in the unshocked xenon is not difiiisive and we rely on the flux-limited diffusion to recover the optically thin 
free-streaming limit. This free-streaming approximation is accurate enough as long as the radiation transport from 
the precursor back to the shock is negligible, in contrast with the almost omnidirectional photon distribution function 
as assumed in the diffusive limit. A fraction of the upstream radiation expands sideways and heats the polyimide 
wall ahead of the primary shock. This will ablate the polyimide of the wall. The resulting inward polyimide flow is 
visible in the Y velocity panel of Fig. |2] This infow extends at 13 ns to x 2100 /im and does have a magnitude of 
about lOkm/s. The exact magnitude of this flow might depend on the radiation transport fidelity used in simulations. 
The ablated polyimide compresses the xenon as can be seen in the density panel in Fig. |2]by the faint tilted feature 
between x » 1700 /vm and x =s 2100 yum. The resulting wall shock has an angle with the primary shock and their shock 
properties were analyzed in 

m. 

The material identity in the top left panel of Fig. |2]also demonstrates the entrainment of xenon in between the 
beryllium and polyimide. In [3J the entrained flow was shown to first get shocked by the wall shock and then again 
shocked near the tripple point of the wall shock and the primary shock. We also mention that in our simulations the 
entrained shear flows produce Kelvin-Helmholtz roll-ups at for instance x « 1400 fim. 

The simulation presented in this section was performed in 3D Cartesian geometry. The axi-symmetry in the 
problem does however allow to perform this simulation in 2D cylindrical rz-geometry as well. We used the same 
settings for the numerical radiation-hydrodynamics solvers. The computational domain is -150 < z < 3900 and 
< r < 900 in microns, where r is the radial coordinate and z is now the coordinate along the tube. The effective 
resolution is 2560 x 512 using two levels of refinement, so that the cell sizes in the rz-geometry correspond to the cell 
sizes in the xy-plane in the 3D Cartesian simulation. The 2D mesh is decomposed in 4 x 4 grid blocks. The mesh 
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Figure 6: Synthetic radiograpli image at 13 ns for the circular nozzle simulation. The X-ray source is located at ix,y, z) = (2,-12, 0) in mm. 



refinement criteria are also the same as for the 3D case. In Fig. |5]the density is shown at the final time 13 ns. The 
results are quite similar to those of the previous simulation. The main difference is in the entrained xenon, which 
requires probably higher resolution to be fully resolved. 

The aim of the CRASH project is to predict certain properties of the compound shock in our high-energy-density 
laser experiments. From the experiment we obtain X-ray radiograph images using He-ff emission from a backlit 
pinhole source transmitting through the experimental target. These images show fundamentally where the dense 
xenon is. To be able to make a direct comparison between the observations and simulations, we added the capability 
to the CRASH code to create synthetic radiograph images. The details of the implementation and verification are 
presented in Appendix B The radiograph for the circular nozzle simulation is shown in Fig. |6] We locate the X-ray 
source at (x, y, z) = (2000, - 1 2000, 0) fim. We have blurred the image to account for the finite pinhole size and finite 
exposure time in the experiment. We also added Poisson noise to mimic the finite photon count. The dense xenon 
behind the primary shock at x = 1700 jum and the wall shock is clearly visible. We can also see the entrained xenon 
and the Kelvin-Helmholtz roll-up to the left of the primary shock. 



3. Elliptical nozzle 

The setup of the elliptical nozzle is very similar to that of the circular nozzle. The main difl'erence is that the wide 
tube of 1200 //m inner radius changes the cross-section down the tube into an ellipse with a major axis of 1200 jum 
and a minor axis of 600 yum. The first 1.1 ns of the simulation is, as for the circular nozzle, performed with the H2D 
radiation-hydrodynamics code to determine the laser energy deposition. For the tube geometry in H2D a straight tube 
is used with a diameter of 1200jt/m. The output of H2D is used to initialize the CRASH code following the recipe 



outlined in Appendix A If we use the notation that x is the coordinate along the tube and y and z are the two directions 
transverse to the tube, then we remap the straight tube of H2D to CRASH using the coordinate transformations shown 
by the equations I A.4I and I A.5 l in which xq - 500 /vm, xi = 750/im, e, =1/2 and s- - 1. That means that the circular 
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Figure 7: The materials in the 3D elliptical nozzle simulation at 1.1 ns (left panel) and 13 ns (right panel). One quarter of the nozzle is shown. The 
color code for the materials is the same as in^ The black iso-surface indicate the primary shock. 

tube is tapered into an elliptical shaft between xq - 500 //m and xi = ISOj-im. The domain size, effective resolution, 
refinement criteria, boundary conditions, and the used numerical scheme are the same as for the circular nozzle. 

The simulation is performed from 1 . 1 ns to 1 3 ns physical time. The number of finite volume cells increased from 
29.6 million initially to about 45.3 million at the end of the simulation. The computational time was 3.5 days on 1000 
processor cores of the HERA supercomputer. The 3D material identification is shown in Fig. [T] The left panel is 
for 1 . 1 ns and the right panel for 1 3 ns. The color code is blue for beryllium, green for polyimide, red for acrylic and 
yellow for gold. Xenon is for convenience not colored in these panels so that the elliptical shaft is visible. At 13 ns 
the beryllium has moved into this shaft and drives a shock in the xenon like a piston. The shock front is indicated 
with a black surface. The edge of the volume occupied by the xenon is colored in white to visuaUze the entrainment 
of xenon between the polyimide and beyllium. 

The left panels of Fig. [8] show the material location, mass density and electron temperature, respectively, in the 
xz-plane at 13 ns. The basic ingredients of the compound radiative shock structure are seen in these panels. The 
primary shock near x ^ 1700 fim in this plane is curved, since the diameter of the laser spot of 800 fim is smaller than 
the major axis, 1200yL(m, of the elliptical shaft. The ripple in the compressed xenon region behind the primary shock 
was analyzed for similar experimental conditions 1 10]. We also find again the tilted wall shock in front of the primary 
shock. The electron temperature panel shows the temperature peak in the equilibration zone behind the primary shock 
and a decreasing temperature behind that in the radiative cooling zone. The right panels of Fig. |8]show the same 
physical quantities in the vertical yz-plane. The shock structure in this plane is quite similar to the results found for 
the circular nozzle. 

The difference in the shock structure and compressed xenon region behind the shock should also be visible in 
the synthetic radiographs. Indeed, in Fig. |9]we show in the left panel the image produced by an X-ray source at 
(x, y, z) - (2000, - 12000, 0) fim directed at the jcz-plane and in the right panel the image produced by an X-ray source 
at (x, y, z) = (2000, 0, 12000) jim directed at the xy-plane. The compressed xenon behind the primary and wall shock 
are found as dark features in these images. Note that the image in the right panel is darker than the image in the 
left panel, since the rays are going through twice as much xenon (the major axis is twice the minor axis). Also note 
that the rippled structure of the compressed xenon layer behind the primary shock is somewhat smoothed out in these 
images. 
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Figure 8: The output of the 3D elliptical nozzle simulation at 13 ns. The state variables indicated in the plot titles are shown in color as a function 
of the (x, z) coordinates given in microns in the y = plane on the left and (x, y) coordinates in the z = plane on the right. The color bars provide 
the scales. The color code of the material levels in the top panels indicates the same materials as in Fig. |2] The black lines in the top panels are for 
the resolution changes, while the black lines in the middle panels indicate the change in material indentity. 
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4. Conclusions 

In this paper we have discussed the simulations of radiative shocks in nozzle shaped tubes. We recover basic 
properties found in laser-driven shock experiments in high-energy-density facilities and from radiation hydrodynamic 
theories: A compressed layer shows up directly behind the primary shock in which the electrons are heated in an 
equilibration zone by the shock wave heated ions. This is followed by a radiative cooling layer. The emitted photons 
propagate upstream and preheat the precursor ahead of the primary shock. The lateral expansion of this radiation 
ablates the plastic nozzle which results in a wall shock. The simulations are shown to be in agreement with simple 
back-of-the-envelope estimates. 

We have used the two-dimensional version of the Hyades code to determine the laser energy deposition during 
the first 1.1 ns of the experiment. This code simulates on a Lagrangian mesh the radiation hydrodynamics based 
on a multi-group flux-limited diff'usion model and flux-limited electron thermal heat conduction. We are currently 
constructing a laser package in the Eulerian CRASH code so that the laser energy deposition and the produced radiative 
shocks in the xenon at later times are self-consistently evaluated with the same code. 

We plan to perform laser-driven shock experiments in three-dimensional circular and elliptical nozzles. The 
comparison of the elliptical nozzle simulations with such experiments will stress test the performance of the code. 
For the validation we can compare the properties of the shock structure and compressed xenon layers from both 
the simulations and experiments. These 3D properties can in the experiments be extracted from dual, orthogonal 
radiography and then compared to the two orthogonal synthetic radiographs from the simulations. 
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Appendix A. Initializing CRASH with Hyades 

The laser energy deposition during the first 1.1 ns is evaluated with the Lagrangian radiation hydrodynamics code 
Hyades 2D (H2D) [4]. The output of H2D is on a distorted logically Cartesian mesh of n^- x n,, cells. When we 
initialize the CRASH code with H2D, this mesh is first triangulated by splitting each quadrilateral along the shorter 
diagonal. The interpolation from the triangulated Hyades grid requires finding the triangle that surrounds the center 
of a given grid cell of the CRASH code. A simple linear search can become very inefficient when we have many grid 
cells (order of a hundred thousand) per processor To accelerate this method, first we create a uniform grid with about 
200 by 200 resolution that covers the whole domain. For each rectangular cell in the uniform grid we find and store 
the list of triangles that intersect it. This can be done very fast. Then we perform the interpolation onto the CRASH 
grid by first finding the rectangular cell that surrounds it, and then we only check the triangles that intersect this cell 
to see which one contains the CRASH grid cell center. 

The H2D simulations are performed on a Lagrangian grid where all cells correspond to a unique material. These 
materials are identified by a material index in the H2D output. Our code uses an Eulerian grid and we track the material 
by means of level set functions. These level set functions are smooth and signed distance functions, initialized with 
the following algorithm. For each cell / of the H2D grid having material index o-,, the level set function for material 
ffj is set to the distance to the closest H2D cell that contains a material different from o-,: 

daiii) - + min |r/ -r;|, (A.l) 

where r, is the location of cell /. The level set functions for the other materials, j3 + a-,, are set to the negative distance 
to the closest cell containing material yS: 

c/sO') = - min |r,- -r;|. (A.2) 
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We interpolate these level set functions to the cells on the CRASH grid. In our first-order level set scheme the material 
in the cell is indicated by the largest level set function. For the simulations with H2D of our shock tube experiments 
we use xenon, beryllium, polyimide, gold, acrylic and "vacuum". In CRASH, vacuum is reassigned as polyimide 
with low mass density. At later times during the evolution, the location of material m follows the simple advection 
equation for the level set function 

+ V ■ {d„,») = d„,V ■ u. (A.3) 

ot 

Also here, the material having the largest d,„ is assigned to be the material of the cell. 

A 3D nozzle is created from 2D Hyades as follows. We first read the H2D shock tube output and triangulate the 
data. After that we transform the coordinates in the transverse directions y and z- The circular cross section of the 
tube shrinks into an ellipse or smaller circle by means of 



y ^ y 



X — X{) 

1 + (ey - l)max 0, min(l, ) 

' Xl- Xq 

, X — Xq 

1 + (e. - l)max|0, mm(l, ) 

Xl - X() 



(A.4) 
(A.5) 



where £, and s, are the factors with which the y and z coordinates contracts for x > xi along the tube. This shrinking 
varies linearly between xq and xi. For x < xq the tube is not modified. This means that the plasma at the far end 
is relocated closer towards the axis. This will change the output of H2D to an unphysical state (the solution of the 
straight tube is not the same as the solution of the nozzle). The impact of this change in geometry is however minimal 
since at 1 . 1 ns the shock dynamics did not yet reach the far end of the shock tube. 



Appendix B. Synthetic radiographs 

The CRASH code has the capability of generating synthetic radiographs during the simulations. The time fre- 
quency of the plotting, the location (or possible locations) of the X-ray source(s) and the orientation, size, and number 
of pixels of the radiograph images are all input parameters. 

In the CRASH shock tube experiments we typically produce radiographs by transmitting 5.18 keV (mainly V He- 
alpha line) X-rays through material whose temperature is of the order of 50 eV. These X-rays interact only weakly 
with the free or outer electrons. We can therefore use cold opacities for the materials in the experiment to evaluate 
the transmission. X-ray photons at this energy experience negligible refraction in the experimental target. So it is 
legitimate for our radiograph calculation to use straight-line analysis. 

These images are line-of-sight plots that calculate the optical depth along the straight lines (rays), connecting the 
X-ray source to the center of the image pixel. The optical depth is an integral of the mass density multiplied by the 
specific opacity characteristic for the material and the spectrum of the X-ray source along the ray: 

D= p(r{l))K(r(r))dh (B.l) 
Jo 

where k is the specific opacity of the materials used in the experiments and the ray is parameterized by the distance / 
to the source. In our shock-tube experiments fs'], we assume that the absorption is dominated by the X-ray line near 
5.18 ke V and the specific opacities of the commonly used materials are 79.4, 0.36 and 2.24 m^/kg for xenon, beryllium 
and polyimide, respectively. In our experiments we also use acrylic and gold, but these materials are outside the view 
of the radiograph, so that we do not need their specific opacity values accurately. For acrylic we use the same value 
as for polyimide and we use the large value 10'" m^/kg for gold. 

The parallel algorithm is implemented in the following way: For each ray and each grid block we determine first 
the segment of the ray if it intersects the block. Then we integrate for each ray and block the optical depth iBTI using 
a trapezoidal rule and a tri-linear interpolation. The step size of the integration is proportional to the cell size of the 
block. Once the integration is done for all blocks, we add up for each ray all integrals over the block segments using a 
call to MPLreduce. For the rz-geometry we consider the 2D blocks as rings in 3D with rectangular cross-section. The 
integration along the ray segments is performed in 3D, but we use a bi-linear interpolation to obtain pK at the required 
integration points. 
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We have verified the implementation with an analytical problem where we integrate along rays through a 3D 
sphere of a given density profile and specific opacity of value one. The density inside the sphere of radius R is defined 
as p = - if the radial distance from the center of the sphere is r < R, while p - for r > R. The density can be 
integrated analytically along a ray passing through this sphere at a distance 6 from the center The total length of the 
segment of this ray inside the sphere is 25 -2 V/?^ - 6^. With these definitions, the optical depth |BT| along this line 
segment becomes 



p(s)ds = (R^-S^- s^)ds = -(R^ 
CO J-s 3 



• 50 ' , (B.2) 



where s is the coordinate along the ray such that \s\ - Vr^ - 6^ is the distance to the middle of the line segment. 
The radiograph image is always a plane through the origin. A pixel with coordinates (x',y') of this image is then 
at a distance d - -\/x'^ + y'^ from the origin (here and in the following the prime indicates the coordinate system 
associated with the radiograph image). The location of the X-ray source is assumed to be on the z'-axis at a distance L 
from the center of the sphere. The ray connecting this source and the pixel {x',y') will then have a minimum distance 
to the center of the sphere of 6 = dL/ yjlT+cP, as can be deduced from similar triangles. The final formula for the 
optical depth along this ray now becomes 



D — max 



»4 



R^- 



1 +L'^l(x'^+y'^) 



3/2^^ 

(B.3) 



In the verification test we place the X-ray source at {x,y,z) - (60,80, 100) and the radiograph image is always 
orthogonal to the fine connecting the source and center of the sphere, so that the distance between the source and 
image is L = 100^/2. The density sphere has a radius R - 10. The left panel in Fig. IB. 101 shows the simulated 
image for a grid resolution of 40^ cells at the base level and one level of refinement for y > 0, z > 0. The right panel 
shows the difference relative to the analytical formula lB.31 Note that the error is smaller where the rays go through 
the refined region of the grid. The grid convergence is seen in Fig. IB.lll for three different resolution at the base level: 
20^*, 40^, and 80^*. The relative error is calculated from 

_ Zii \Di - Aef;l 
" " V/ in I ' ^^-^^ 



where D is the simulated optical depth, D,ef is the analytical reference IB . 3 1 and / = 1, . . . / indexes the pixels of the 
image. We obtain second-order convergence. We have also verified the implementation for the rz-geometry, in which 
case the sphere is a circle rotated around the symmetry axis. 

The simulated radiographs can be compared with experimental backlit pinhole radiographs for validation studies. 
To make these images appear more similar, we need to take into account the finite pinhole size of about 20 mm 
diameter, the finite exposure time of about 0.2 ns and the effect of a finite number of collected photons (typically 50 
photons per 100 nm- image pixel). The first two effects are approximated by smoothing the synthetic radiograph over 
a few pixels. The finite photon count can be taken into account by using 

D' = - ln[P(50)e"^] = D - lnP(50), (B.5) 

instead of D, where P(50) is a random number with a Poisson distribution and a mean value of 50. The second 
equality shows that this is the same as subtracting the logarithm of these random numbers from the original simulated 
radiograph during the post-processing step. In Fig. IB. 121 the original synthetic radiograph is compared with the 
post-processed image. 



Appendix C. Improved diffusion operator at resolution changes 

In this appendix we present an improved conservative and spatially second-order implicit scheme for the radiation 
diffusion and heat conduction. For convenience, we only derive this scheme for the heat conduction. The generaliza- 
tion for radiation diffusion is straightforward. 
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Figure B.IO: The radiograph image for an analytical density sphere on a non-uniform mesh and a specific opacity with value one. The X-ray 
source is located at (x,y, z) = (60, 80, 100). The left panel shows the simulated image, while the right panel is for the error relative to the analytical 
solution. 
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Figure B.ll: The LI error of the optical depth relative to the analytical reference for 20^, 40^ and 80^ base resolutions and one level of refinement 
for .V > 0, >■ > and z > 0. 
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Figure B.12: Original syntlietic radiograph (left panel) and the processed image after smoothing and adding Poisson noise to mimic the experimental 
radiograph (right panel). 



Discretizing the electron thermal heat conduction implicitly in time leads to the linearized backward Euler equation 
for the electron temperature T 

c*, — - — = V ■ c:vr+', (c.i) 

where Cve is the electron specific heat and Ce is the heat conduction coeifcient. The time level * corresponds to the 
state before the implicit update. During the implicit advance with time step Af, the coefficients are frozen in at time 
level *, resulting in a temporally first order scheme in general. This equation can be recast in a linearized equation for 
the change Ar = r"+' -T*: 

C* 



AT = v -cyr. (C.2) 



Af 

The right-hand-side depends only on time level *. Once Eq. IC.2l is solved using a Hnear solver, the electron energy 
density can be updated using 

= e; + c^,(r'+' - T*), (C.3) 

to conserve the energy. 

A discrete set of equations is obtained by applying a finite volume method to Eq. IC.2I To make this scheme 
spatially second-order accurate on a uniform mesh, we need a second order accurate thermal heat flux at the face 
centers. This is achieved by approximating the gradient of the electron temperature with a central difference using the 
ceU-centered values 

Ti-Tj 



X 



V ■ (CyT)dV = VSijCij-^ (C.4) 



where the control volumes are indexed by /, each having a volume V,-. The index j is for the neighboring cells having 
a common interface with area S ij and a distance between the cell centers is |x, - Xj\. The heat conduction coefficient 
at the face is the arithmatic average of the coefficient at the two neighboring cell centers, Ceij = (C^/ + Cej)/2. 

To obtain a second-order heat flux at the resolution change, we need a third-order interpolation of the temperature 
in the ghost cells. Such an interpolation was previously used in the context of Hall magnetohydrodynamics (MHD) 
I12I1 . This interpolation is only needed for the fine cells, since the flux at the coarse side will be obtained as the sum 
of the fluxes at the neighboring fine cells to preserve conservation of the scheme (]3\. The implementation for the 
heat conduction is different from the Hall MHD, since for the heat conduction we also have to maintain positivity 
of the temperature and avoid spurious oscillations. For convenience we will restrict the analysis to two-dimensional 
domains. The third-order interpolation for the temperature value at the fine ghost cell (0, j), indicated by the dashed 
circle, is first performed along the coarse cell values in the transverse direction to obtain the temperature at (-1/2, j) 
as depicted in Fig. IC.13I 

5r_i/2j-3/2 + 30r_i/2j+i/2 - 3r-I/2,;+5/2 

^-i/2j = • (<--5) 
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Figure C.13: The coefBcients of the third-order interpolation near the resolution changes (similar to fl^)- The order of the interpolation is first 
along the coarse grid cells tangential to the resolution change, followed by an interpolation along the fine cells to obtain a third-order temperature 
value in the fine ghost cell indicated by the dashed circle. 



To guarantee the positivity of the interpolated temperature, the value is clipped by the maximum and minimum values 
of the surrounding points 

r_i/2j = max{min[r_i/2j, max(r_i/2j-3/2, 7'-i/2j+i/2)], (C.6) 

min(r_i/2j-3/2, 7^-1/2,7+ 1/2))- (C.7) 

The value in the fine ghost cell can now be obtained by a parabolic interpolation along the fine cells in the direction 
normal to the refinement interface 

8r_i/2j+iori.,-3r2j 

15 

For this interpolation we need again to clip the obtained value with the surrounding temperatures 



Oi-l/2J-t-iUil.j--?i2j 
■'(),) - — ■ (^-o) 



Toj = max{min[ro,;, max(r_i/2,;, Tij)],mm(T-i/2,j, Tij)}. (C.9) 



This ghost cell value is used in the heat conduction formula lC4l We still need a second order accurate heat conduction 
coefficient Cdj at the fine face center. This amounts to a second-order prolongation of the heat conduction coefficient 
to obtain the ghost-cell values at the fine side, followed by an averaging of the cell centered coefficients to the face 
center Once the thermal heat fluxes at the fine side are obtained, conservation is restored by copying the fine fluxes to 
the coarse side at the resolution changes. Note that this scheme is different from |14] where conservation of the flux 
at the resolution changes is enforced in the strong sense by enforcing the flux on the coarse side to be equal to each of 
the fluxes on the fine side. 



In the analysis we assumed a Cartesian mesh. Generalization to curvilinear grids is presented in II12I1 . The same 
third-order interpolation procedure with clipping is also used for heat conduction along the magnetic field lines in the 
context of solar wind modeling U5|] . 
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Figure C.14: Heat conduction test in rz-geometry with uniform coefficients. The top panel shows the electron temperature in color contour at the 
final time t = 3/2. The resolution changes are indicated by the black box. The middle panel is for the temperature error relative to the reference 
solution, while the bottom panel shows the relative enor of the original heat conduction scheme 1(5(1. 



To verify that the improved implicit heat conduction and radiation solver is second order, we first demonstrate a 
uniform heat conduction test in rz-geometry. This test follows the time evolution of the electron temperature profile 
for a purely heat conductive plasma similar as described in [SJ. We set the electron specific heat to be one and further 
assume the heat conductivity Ce to be a constant. The time evolution for this problem is a product of a Gaussian 
profile in the z-direction and the Bessel function Jq in the r-direction: 

T = r™„ + To ^/^ e-^Jo(br)e-''''^'', (C.IO) 
^J4■nCet 

where b ~ 3.8317, T^i„ = 3, To = 10, and C, = 0.1. 

The domain size for this problem is -3 < z < 3 and < r < 1, which we decompose in 3 x 3 grid blocks of 30 x 30 
cells each. The central block, -1 < z < 1 and 1/3 < r < 2/3, is mesh refined by one level. At the outer boundary 
we fix the ghost cells to the exact solution IC.IOI except for r = 0, in which case we apply a symmetry condition. 
We simulate the time evolution from f=ltof = 3/2 using the GMRES iterative solver with a Schwarz-type ILU 
preconditioner To achieve second-order time integration we use the Crank-Nicolson approach. The latter is possible 
since the coefficients in the problem are all temporally invariant. 

The final solution with 90^ base resolution is shown in the top panel of Fig. IC.14I The color is for the electron 
temperature. The change in the grid resolution with one level of refinement is indicated by the black line. To demon- 
strate that the temperature error at these resolution changes is not significantly larger than at the uniform part of the 
mesh, we also plot the spatial distribution of the temperature error in the bottom panel. The maximum error is clearly 
at the coarse grid away from the resolution changes. 

We use the relative maximum error to study the grid convergence. This error is defined by 

max,=i,..../ IF; - r,-ef,| 
max,=i,...,/ ir^ef/l 

in which T is the numerical solution on the ; = 1 , . . . , / control volumes and Tref is the analytical reference in Eq. 
IC.IOI The second-order convergence rate is depicted in Fig. IC.lSl for grid resolutions of 90^, 180^, 360^, and 720^ at 
the base level. 

To demonstrate that the new heat conduction and radiation solver does not under- or over-shoot near discontinu- 
ities, we show the Mach 5 non-equilibrium gray-diffusion test of jl^. This test uses non-uniform radiation diffusion 
coefficient and Planck opacity that depend on the density and temperature as defined by Dr - 0.0175(yr)^^-/p and 
CKp = 10^ /Dr- In [5], this verification was transformed to a heat conduction test. Here, we will keep it as a radiation 
test. We add a Mach -5 flow to the ID initial condition, so that both the shock and radiation condition will move to the 
left with a Mach -5 velocity. This solution is then rotated counter-clockwise over an angle tan '(1/2) on a 2D grid to 
make the test problem more difficult. The domain size is -0.0384 < x < 0.0384 by -0.0048 <y < 0.0048 and within 
the region -0.0128 < x < 0.0128 and -0.0016 <y < 0.0016 we refine the mesh by one level. The grid is designed in 
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Figure C.15: Relative maximum error for the heat conduction test in r^-geometry and with uniform coelficients. Second-order convergence is 
achieved on a non-uniform grid. 

such a way that during the time evolution both the shock and the radiation precursor front travel through a resolution 
change. For the simulation of the hydrodynamic part, we use the HLLE scheme with the generalized Koren limiter by 
setting jS = 3/2 and a CFL of 0.8. The radiation diffusion and energy exchange between the radiation and material are 
solved implicitly with the new implicit method using the GMRES iterative solver with a block ILU preconditioner 

The material and radiation temperatures at the final time, t = 0.0025, are shown in the left and right panel of Fig. 
IC.16I respectively. The solid, dotted, dashed, and dot-dashed black lines indicate the base resolutions of 192 x 24, 
384 X 48, 768 x 96, and 1536 x 192, respectively. The blue line is the semi-analytic solution of 1 16]. The new scheme 
correctly solves the precursor front at x » -0.022 and shock front at x ^ 0.0085 without spurious oscillations. 
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Figure C.16: Material (left panel) and radiation (right panel) temperature at the final time for the Mach 5 radiative shock tube problem rotated on a 
non-uniform grid. The blue line is the reference solution obtained from |16j. The insets show the temperature jump at the hydro-shock (left panel) 
and the radiation precursor front (right panel). 
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